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The basic equations are derived for compressible flow in a 
labyrinth seal. The flow is assumed to be completely turbulent 
in the circumferential direction where the friction factor is 
determined by the Blasius relation. Linearized zeroth and first- 
order perturbation equations are developed for small motion about 
a centered position by an expansion in the eccentricity ratio. 
The zeroth-order pressure distribution is found by satisfying the 
leakage equation while the circumferential velocity distribution 
is determined by satisfying the momentum equation. The first- 
order equations are solved by a separation of variables solution. 
Integration of the resultant pressure distribution along and 
around the seal defines the reaction force developed by the seal 
and the corresponding dynamic coefficients. The results of this 
analysis are compared to published test results. 


INTRODUCTION 


The problem of self excited vibration in turbomachinery due 
to labyrinth seals has led to the development of many analyses 
which attempt to model the physical phenomenon so that the 
problem can be better understood and therefore solved. The 
shortcoming with the analyses which have been presented to date 
is that they are difficult to understand and require limiting 
assumptions such as ignoring the area derivative in the 
circumferential direction, assuming that the friction factor is 
the same for all surfaces, and assuming that the flow coefficient 
is constant along the seal. These assumptions may be of some use 
mathematically, but do very little for the understanding of the 
physical occurrence. The first steps toward analysis of this 
problem were taken by Alford [1] , who neglected circumferential 
flow and Spurk et al [2] who neglected rotation of the shaft. 
Vance and Murphy [3] extended the Alford analysis by introducing 
a more realistic assumption of choked flow, Kostyuk [4] 
performed the first comprehensive analysis, but failed to include 
the change in area due to eccentricity which is responsible for 
the relationship between cross-coupled forces and parallel rotor 
displacements. Iwatsubo [5,6] refined the Kostyuk model to show 
this relationship by including the time dependency of area 
change, but he neglected the area derivative in the 
circumferential direction. Kurohashi [7] incorporated dependency 
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of the flow coefficient on eccentricity into his analysis, but 
assumed that the circumferential velocity in each cavity was the 
same . 


The analysis presented here includes the variation of the 
area in the circumferential direction due to eccentricity and 
incorporates as many of the physical phenomena in the flow field 
as was thought necessary to produce an adequate result. The main 
purpose of this paper is to present a unified and comprehensive 
derivation of a reduced set of equations and a new solution 
format for those equations. The results of this analysis are 
compared to the published test results of Wachter and Benckert 
[8,9,10]. 


NOMENCLATURE 


A- Cross sectional area of the cavity (L**2 ) ;def ined in text 
B- Height of labyrinth seal strip (L);defined in fig.(l) 

C- Direct damping coefficient (Ft/L) 

Cr- Nominal radial clearance (L) ; defined in fig.(l) 

D- Hydraulic diameter of cavity (L) ; introduced in Eq.(3) 

H- Radial clearance (L) 

K- Direct stiffness coefficient (F/L) 

L- Pitch of seal strips (L); defined in fig.(l) 

NT- Number of seal strips 
P- Pressure (F/L* *2) 

R- Gas constant 

Rs- Radius of seal (L) ; Defined in fig.(l) 

T- Temperature (T) 

Rw- Surface velocity of rotor (L/t) 

V- Average velocity of flow in circumferential direction (L/t) 
a,b- Radial seal displacement components due to elliptical whirl 
(L) ; defined in Eq.(13) 

ar- Dimensionless length upon which shear stress acts on rotor 
as- Dimensionless length upon which shear stress acts on stator 
c- Cross coupled damping coefficient (Ft/L);in Eq.(18) 
k- Cross coupled stiffness coefficient (F/L); in Eq.(18) 
m- Leakage mass flow rate per circumferential length (M/Lt) 
mr, nr, ms, ns- Coefficients for Blasius relation for friction 
factor; defined in Eq.(3) 

t- Time (t) 

C 0- Shaft angular velocity (l/t) 

P~ Density of fluid (M/L**3) 

V- Kinematic viscosity (L**2/t) 

€= e/Cr 

Y- Ratio of specific heats 

Kq - Dimensionless cross-coupled stiffness parameter ; defined in 
# Eqs . ( 27 ) 

Eq- Dimensionless entry-swirl parameter; defined in Eqs. (27) 


Subscripts 

o- Zeroth-order component 
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i- i-th chamber value 
1- First-order component 
x- X-direction 
y- Y-directin 
r- Reservoir value 
s- Sump value 


PROCEDURE 


The analysis presented here is based on the see-through type 
of labyrinth seal shown in figure 1 . The continuity and momentum 
equations will be derived for a single cavity control volume as 
shown in figures 2 , 3 , 4 , and 5. A leakage model will be employed 
to account for the axial leakage. The governing equations will 
be linearized using perturbation analysis for small motion about 
a centered position. The zeroth-order continuity and momentum 
equations will be satisfied to yield the steady state pressure 
and velocity for each cavity. The first-order continuity and 
momentum equations will be reduced to linearly independent, 
algebraic equations by assuming an elliptical orbit of the shaft 
and a resulting harmonic response for the pressure and velocity 
fluctuations. The force on the shaft will be found by 
integration of the pressure fluctuations along and around the 
shaft. Using the equations for forced motion of the shaft, the 
stiffness and damping coefficients will be found. 


ASSUMPTIONS 


1 ) Fluid is considered to be an ideal gas . 

2) Pressure variations within a chamber are small compared to the 
pressure difference across a seal strip. 

3) The frequency of acoustic resonance in the cavity is much 
higher than that of the rotor speed, 

4) Added mass terms are neglected. 

5) The eccentricity of the rotor is small compared to the radial 
seal clearance. 

6) In the determination of the shear stresses in the 
circumferential direction, the axial component of velocity 
is neglected. 

7) The contribution of shear stress to the stiffness and 
damping coefficients is neglected. 

GOVERNING EQUATIONS 


Continuity Equation 


Referring to the control volume in figures 2 and 
continuity equation for the control volume shown is: 
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where the transverse surface area Ai is defined by; 

A ± - [B i + (Cr+eH 1 ) i + B i+1 + (Cr+eHi) i+i ]L ± /2 

Momentum Equation 

The momentum equation (2) is derived using figures 4 and 5 which 
show the pressure forces and shear stresses acting on the 
control volume. This equation includes the area derivative in 
the circumferential direction, which was neglected by Iwatsubo 
[5,6]. 
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where ar and as are the dimensionless length upon which the shear 
stresses act and are defined as; 


as = Lj/Lj ar = ( 2B+Lj )/Lj 

for teeth on the rotor and as; 

as = (2B+Lj )/Lj ar = Lj/Lj 

for teeth on the stator. Blasius [11] determined 
stresses for turbulent flow in a smooth pipe could 
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that the shear 
be written as ; 


(3) 


This relationship is applied to the labyrinth surface for the 
hydraulic diameter; 

D = 2 (Cri+B)Li 
(Cri+B+Li) 


and where Urn is the mean flow velocity relative to the surface 
upon which the shear stress is acting. The constants mo and no 
can be empirically determined for a given surface from pressure 
flow experiments. However, for smooth surfaces the coefficients 
given by Yamada [12] for turbulent flow between annular surfaces 
are : 

mo = -0.25 no= 0.079 

For the control volume shown in figure 4, ns and ms represent 
tne constants for the stator surface and nr and mr represent 
those for the rotor surface. This feature allows one to tailor 

this analysis to a particular surface such as honeycomb or an 
abradable one . 
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Substituting the mean flow velocity relative to each surface 
and referring to figure 5, which accounts for the pressure on the 
sides of the control volume due to the eccentricity of the seal, 
the circumferential momentum equation for the seal is: 
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If Eq.(l) times the circumferential velocity is now subtracted 
from Eq. (4) , the following reduced form of the momentum equation 
is obtained: 
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In order to reduce the number of variables, all of the 
terms are replaced with pressure terms using the ideal 


( 6 ). 


P ± * PRT 


density 
gas law 

( 6 ) 


Furthermore, in order to make the perturbation analysis easier, 
the following substitution is made in the continuity equation; 


Leakage Equation 


• 2 • 2 

in.. - m, = i+r i 

1+1 i — 


To account for the leakage mass flow rate in the continuity 
and momentum equations, the leakage model of Neumann [13] was 
chosen. This model predicts leakage and pressures fairly 
accurately and has a term to account for kinetic energy 
carryover. However, the empirical flow coefficient relations 
given by Neumann were discarded in favor of the equations of 
Chaplygin [14] for flow through an orifice. This was done to 
produce a different flow coefficient for succeeding contractions 
along the seal as has been shown to be the case by. Egli [15]. 
The form of the model is: 
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where the kinetic energy carryover coefficient fJ 2 is defined for 
straight through seals as? 


where 


V2 ‘4 


NT 
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and is unity, 
groove seals. 


by definition, for interlocking and combination 
The flow coefficient is defined as; 

Y-l 
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For choked flow, Fliegner ' s formula will be used for the last 
seal strip. It is of the form; 


where 
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LINEARIZATION 


Since an analytical solution to the governing equations 
not available, the continuity and momentum equations will 
expanded in the following perturbation variables: 


P i = P oi + eP > 
V i * V oi + eV, 


H i " <c ri + eiy 

A i ’ A oi + “A 


is 

be 


where £ =e/Cr is the eccentricity ratio. The zeroth-order 
equations define the leakage mass flow rate and the velocity 
distribution for a centered position. The first-order equations 
define the perturbations in pressure and circumferential velocity 
due to a radial position perturbation of the rotor. Strictly 
speaking, the results are only valid for small motion about a 
centered position. 

Zeroth-Order Solution 

Eq. (9) is used to determine the pressure distribution along 
the shaft in the following manner. The leakage is solved using 
Eq.(7) or Eq.(8), depending on the operating conditions. To 
determine if the flow is choked or not, assume that the pressure 
in the last cavity is equal to the critical pressure for choking. 
Using this pressure, find the leakage from Eq. (8) and then use 
Eq.(7) to determine the reservoir pressure necessary to produce 
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this condition. Based on this pressure, a determination can be 
made whether the flow is choked or not. The associated pressure 
distribution is determined by employing the correct leakage, 
along with a known boundary pressure, and solving Eq. (7) one 
cavity at a time. 
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(9) 


For cavity i, the zeroth-order circumferential momentum equation 
is: 
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With the pressure determined, the only variables remaining in the 
momentum equation for the cavities are the velocities. Given an 
inlet tangential velocity, a Newton root finding approach can be 
taken whereby Eq.(10) is solved for the i-th velocity, one cavity 
at a time. This is done starting at the first cavity and working 
down stream. 

First-Order Solution 

The governing first-order equations (11,12), define the 
.pressure and velocity fluctuations resulting from the seal 
clearance function. The continuity equation (11) and momentum 
equation (12) follow: 
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where the X's and G's are defined in Appendix A. If the shaft 
center moves in an elliptical orbit, then the seal clearance 
function can be defined as: 


eHj a -a CO suit cos0 -b sinwt sinQ 


(13) 



The pressure and velocity fluctuations can now be stated in the 
associated solution format: 


*4* •>!«• 
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4- 4 . 

v 1 i = V ci c °s (0+wt) + V gi sin(0+O)t) + V~,.cos(0-u)t) + V gi sin(0^t) (15) 

Substituting Eqs. ( 13 ) , ( 14) , and (15) into Eqs.(ll) and (12) and 
grouping like terms of sines and cosines (as shown in Appendix B) 
eliminates the time and angular dependency and yields eight 
linear algebraic equations per cavity. The resulting system of 
equations for the i-th cavity is of the form; 
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The A matrices and column vectors B and C are given in Appendix 
B. To use Eq.(16) for the entire seal solution, a system matrix 
must be formed which is block tridiagonal in the A matrices. The 
size of this resultant matrix is (8(NT-1)) x (8(NT-1)) since 
pressure and velocity perturbations at the inlet and the exit are 
assumed to be zero. This system is easily solved by various 
linear equation algorithms, and yields a solution of the form: 
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DETERMINATION OF DYNAMIC COEFFICIENTS 


The force-motion equations for a labyrinth seal are assumed 
to be of the form: 
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The solution of Eq.(18) for the stiffness and damping 
coefficients is the. objective of the current analysis. For the 
assumed elliptical orbit of Eg. (.13)* the X and Y components of 
displacement and velocity are defined as : 


X = a coswt 
Y * b sinwt 


X - -aw sinwt 
Y = bu) coswt 


Substituting these relations into (18) yields: 

-3? = -Ka coswt - kb sinwt + Caw sinwt - cbw coswt 

X (1< 
-F = -ka coswt - Kb sinwt - caw sinwt - Cbw coswt 

Redefining the forces/ Fx and Fy, as the following? 

F = F coswt + F sinwt 

x XC XS / 2( 

F » F coswt + F sinwt ^ 

y yc ys 

and substituting back into (19) yields the following relations: 


-Fxc ■ Ka + cbw 
— Fyc ®-ka + Cbw 


-Fxs ®-CaW+kb 
-Fys ® Kb + caw 


The X and Y components of force can be found by integrating the 
pressure around the seal as follows: 
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Only one of these components needs to be expanded in order to 
determine the dynamic coefficients. For this analysis, the X 
component was chosen. Substituting Eq. (14) into (22) and 
integrating yields: 
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Substituting from Eq.(17) and (19) into Eq. ( 24) and equating 
coefficients of sinmt and cos cot yields: 


NT-1 , _ -1* - m 
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Equating the definitions for Fxs and Fxc provided by Eqs. (21 ) and 
(25) and grouping like terms of the linearly independent 
coefficients a and b yields the final solutions to the stiffness 
and damping coefficients: 
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SOLUTION PROCEDURE SUMMARY 


In review/ the solution procedure uses the following 
sequential steps; 

a) Leakage is determined from Eq. (7) or (8). 

b) Pressure distribution is found using Eq. (7) . 

c) Velocity distribution is determined using Eq. ( 10) . 

d) A system equation is formed and solved using the cavity 
Eq. (16) . 

e) Results of this solution as defined in Eqs. (17), are 
inserted into Eq.(26) 


RESULTS 


To compare the present analytic solution with the 
experimental results of Wachter and Benckert [8,9,10], the 
following dimensionless parameters are introduced. The 
dimensionless cross-coupled stiffness and entry-swirl parameters 
are defined by Wachter and Benckert as: 
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All of the results presented for comparison in this paper are for 
a seal ■with teeth on the stator, with entry swirl and no shaft 
rotation. Although Wachter and Benkert published results for 
shaft rotation, the data for the operating conditions and seal 
geometry were insufficient for use in this study. The results in 
figures 7, 8, and 9 are from [8] and show the relationship between 
cross coupled stiffness and the entry swirl, for a seal with 
strips on the stator and the geometry shown in figure 6. The 
line shown is the experimental result and the symbols are the 
results from this analytical model. These figures show that the 
model compares favorably to the experimental results in magnitude 
and the overall trend for various operating conditions. The 
figures also show that the model does not yield a consistently 
high or low result. Instead, the model tends to over predict the 
value of the stiffness for a large number of strips and under 
predict stiffness for a small number of seals. This trend is 
probably due to errors in calculating the zeroth-order pressure 
distribution using the leakage model. 

The results in table (1) are from [9,10] for a seal with 
strips on the stator. The results show the effect of change in 
seal parameters such as pitch, number of teeth, radius, strip 
height, and clearance on the cross coupled stiffness. The model 
accurately shows the increase in cross-coupled stiffness due to 
decrease in clearance and decrease in strip height, but it fails 
to remain constant for change of pitch and consistantly over 
estimates the cross-coupled stiffness for the larger radius cases 
by about 26%. 


CONCLUSION 


A clear and understandable analysis utilizing reduced 
equations has been presented for the problem of calculating 
rotordynamic coefficients for labyrinth seals. This paper was 
developed to provide a less restrictive analysis and a’ better 
explanation of the current analyses. The model developed gives 
results that were within 25% of the experimental results 
available. However, this error must be balanced against the 
known uncertainties in the experimental data. This is especially 
important since all of the data used are for a nonrotating shaft 
and the only influence on the cross coupled stiffness was the 
entry swirl. Although Wachter and Benckert published data for a 
rotating shaft, the data were not sufficient to calculate a 
result. Also, the only data available for see-through labyrinths 
is for the type with strips on the stator. For a more rigorous 
test of this and other models, more complete data are required 
oyer a wider range of parameters for different seal geometries. 
Finally, this analysis is only considered valid for the see- 
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through type of labyrinth seal since the model fared very poorly 
in comparison with interlocking and grooved seal data. 
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TABLE I. - COMPARISON OF DATA FOR VARIOUS' GEOMETRIES AND OPERATING CONDITIONS 
FOR SEAL WITH TEETH ON STATOR AND NO SHAFT ROTATION [9,10] 
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Figure 4. - Forces on control volume. 



Figure 5. - Forces on control volume. 


271 




Figure 6. — Configuration used for experiment. 
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APPENDIX A 


DEFINITION OF FIRST-ORDER CONTINUITY 
AND MOMENTUM EQUATION COEFFICIENTS 


G, = 


oi 

RT 


P .L 

r = 01 

h 2 RT 


2 2 

hJl 

. RT 
m 

o 


. y, (4s - 5) 
J oi (<&i * °i) - ^ — T 


oi 


M 


G 4 " 


2 2 

- y i U 2 „ 2 


m RT 
o 


p 0l-l + ™ ^ <*-« 


L P oi 


2 2 

G,- = • ■ Cc - P ; , , - + — y, (4s-5) 

5 m RT x+1 oi+1 7T 

o 


Aoi 

-V; 

y oi + 


h-A 

( p ° i - i \ 

\ Y / 


i h 

■ - 1\ / 

'oi+l V 

Y ) { 


Y - 1 
Y 


oi - 1 
P oi 


~y 

Y 


Y - 1 
Y 


G, = m_ 


°i - °S+r 


6 


P .A . 
v _ oi oi 

X 1 " RT 


T asL(2+ms) t. arL(2-hnr) 


X 0 = m H — — 

2 o Vox 


(Rw-Voi) 


_ vf; + , v v . p + mo ,, 

*3 P * P . RTtho V oi -V oi-l P oi + ir y i ( 4s-5) 
ox oi 


1 /Y - V 


P oiV Y 


oi-1 

*oi 


Y - 1 
Y 


(V oi - W 
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— (V • 
TT v oi 


- V oi-l> (4s * 5) 
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APPENDIX B 


SEPARATION OF CONTINUITY AND MOMENTUM EQUATIONS 
AND DEFINITION OF SYSTEM MATRIX ELEMENTS 


CONTINUITY: 


cos ( 0 +tot) 


G . P + . (to+^p-) + G V + G„ P + + G, P + 

1 sx y Rs J 1 Rs sx 3 ci 4 ci-1 


(a-b) 


sin (e-hot) : -G, Pt, + toj- Gn — V+ + G 0 P^ + 0 A P? 4 _, + G„ P + 


G 6 / V oi 
2 \ \ 


1 ‘ci y Rs ■ ur r "1 Rs ci ' si u 4 si— 1 ’ “5 si+1 
+to) (b-a) 


cos (0-tot) : G 


1 \ Rs 


Voi to) P" + G, V . + G 0 P . + G, P . . +G, P . 

si 1 Rs si 3 ci 4 ci-1 5 ci+I 


(a+b) 


sin (9-OJt) : G, («. - P^ - + Gj P^ + G 3 P^ + G, P^ 

+ G e P_, - -4 (u (a+b) 


•5 si+1 2 V Rs 


MOMENTUM: 


Voi) „+ x ji B + 




cos ce+0)t); X, (a. + V al + P sl + X 2 V ci -mo + X 3 


+ x 4 F ci-i ■ r < b - a) 


\ A 

sin (6+ct) : - Xl (« + ^i) V+ ± - -2i P+ ± + X 2 -mo + X 3 P+ ± + X 4 P^ - 0 


cos (0-Hot) : X, (~i - to) v" + “• P ^ + X„ v", -mo V". , + X„ P“ . = -^ 

1 \ Rs / si Rs si 2 ci cx-1 3 ci-1 2 


(a+b) 


sin (0-tot) : X., ( to - ~i) P 


n , + . 21 p“ + x„ V~ - riio V" . + X, P . + X. P". . = 0 

Rs / ci Rs ci 2 si si-1 3 sx 4 si-1 
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A i _ 1 MATRIX 


a l,2 = a 2,l = a 3,4 = a 4,3 = G 4 

a 5,2 = a 6,l = a 7,4 = a 8,3 = X 4 

a 5, 6 = a 6,5 = a 7, 8 = a 8, 7 = 

The remaining elements are zero 


A i MATRIX 


a l,l G 1 ( W + Rs 


'V 

a„ „ = G, / oi 
3,3 1 to 


a 2,2 “ 

“ G 1 

(- * £ 

A 4,4 = 

1 

o 

H* 

ft- 


a l,2 * a 2, 1 ~ a 3,4 " a 4,3 G 3 


a 5,2 a 6.1 ~ a 7,4 “ a 8.3 X 3 


oi 

5,1 “ a 7,3 " Rs 


a„ „ = a 


-A . 
oi 


6,2 8,4 Rs 


“5,5 - X 1 + 


__ | , oi 

a 6,6 " - X 1 l u + ‘l£ 


V 


a, , = X 


oi 


7,7 1 \ Rs 


- to 


v [ oi 
a 8, 8 X 1 1“ _ Rs 


a 5,6 “ a 6, 5 a 7,8 a 8, 7 X 2 


a l,5 a 3,7 = G 1 Rs 
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a 2,6 a 4,8 " G 1 Rs 

The remaining- elements are zero 


A i+1 MATRIX 


a l,2 ' a 2,l ‘ a 3,4 ■ a 4,3 ■ G 5 


THE REMAIHING ELEMEHTS ARE ZERO 


B and C column Vectors 



C = 





0 


279 



